library(ggplot2)
install.packages("mixtools")
library(mixtools)
library("readxl")
library(lme4)
install.packages("glmm")
library(glmm)
install.packages("gee")
library(gee)

setwd("") # <---- insert your own path to the cytokine, MMP, and cartilage data file here
Lum_data <- read.csv("cartilage_cytokine_MMP_raw_data_for_PLOS_ONE.csv", header = TRUE)
Lum_data <- Lum_data[complete.cases(Lum_data),]
#Which values are greater than 0.0
replace_with_one <- which(Lum_data[,c(4:22)] > 0.0, arr.ind = TRUE)
replace_with_one[,2] <- replace_with_one[,2] + 3
Lum_data[replace_with_one] <- 1.0
Lum_data[Lum_data==-1] <- 0

# Good: [1,5]
# Fair: [6,8]
# Bad: [9,Inf)
Good_rows <- Lum_data[which(Lum_data$Group == "Good"),c(4:22)]
par(mfrow=c(3,4))
hist(Good_rows$GM.CSF..14.,main = "GM-CSF",xlab="Cytokine Presence",col = "green")
hist(Good_rows$IL.1.alpha..22.,main = "IL-1a",xlab="Cytokine Presence",col = "green")
hist(Good_rows$IL.RA..29.,main = "IL-Ra",xlab="Cytokine Presence",col = "green")
hist(Good_rows$IL.2..35.,main = "IL-2",xlab="Cytokine Presence",col = "green")
hist(Good_rows$IL.4..39.,main = "IL-4",xlab="Cytokine Presence",col = "green")
hist(Good_rows$IL.6..43.,main = "IL-6",xlab="Cytokine Presence",col = "green")
hist(Good_rows$IL.8..46.,main = "IL-8",xlab="Cytokine Presence",col = "green")
hist(Good_rows$IL.10..52.,main = "IL-10",xlab="Cytokine Presence",col = "green")
hist(Good_rows$IL.12..57.,main = "IL-12",xlab="Cytokine Presence",col = "green")
hist(Good_rows$IL.18..61.,main = "IL-18",xlab="Cytokine Presence",col = "green")
hist(Good_rows$TNF.alpha..76.,main = "TNF-a",xlab="Cytokine Presence",col = "green")

par(mfrow=c(2,4))
hist(Good_rows$MMP.1..18.,main = "MMP-1",xlab="Cytokine Presence",col = "green")
hist(Good_rows$MMP.2..21.,main = "MMP-2",xlab="Cytokine Presence",col = "green")
hist(Good_rows$MMP.3..27.,main = "MMP-3",xlab="Cytokine Presence",col = "green")
hist(Good_rows$MMP.7..30.,main = "MMP-7",xlab="Cytokine Presence",col = "green")
hist(Good_rows$MMP.9..33.,main = "MMP-9",xlab="Cytokine Presence",col = "green")
hist(Good_rows$MMP.10..37.,main = "MMP-10",xlab="Cytokine Presence",col = "green")
hist(Good_rows$MMP.12..44.,main = "MMP-12",xlab="Cytokine Presence",col = "green")
hist(Good_rows$MMP.13..48.,main = "MMP-13",xlab="Cytokine Presence",col = "green")


Fair_rows <- Lum_data[which(Lum_data$Group == "Fair"),c(4:22)]
par(mfrow=c(3,4))
hist(Fair_rows$GM.CSF..14.,main = "GM-CSF",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$IL.1.alpha..22.,main = "IL-1a",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$IL.RA..29.,main = "IL-Ra",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$IL.2..35.,main = "IL-2",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$IL.4..39.,main = "IL-4",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$IL.6..43.,main = "IL-6",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$IL.8..46.,main = "IL-8",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$IL.10..52.,main = "IL-10",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$IL.12..57.,main = "IL-12",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$IL.18..61.,main = "IL-18",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$TNF.alpha..76.,main = "TNF-a",xlab="Cytokine Presence",col = "gray")

par(mfrow=c(2,4))
hist(Fair_rows$MMP.1..18.,main = "MMP-1",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$MMP.2..21.,main = "MMP-2",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$MMP.3..27.,main = "MMP-3",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$MMP.7..30.,main = "MMP-7",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$MMP.9..33.,main = "MMP-9",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$MMP.10..37.,main = "MMP-10",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$MMP.12..44.,main = "MMP-12",xlab="Cytokine Presence",col = "gray")
hist(Fair_rows$MMP.13..48.,main = "MMP-13",xlab="Cytokine Presence",col = "gray")


Bad_rows <- Lum_data[which(Lum_data$Group == "Bad"),c(4:22)]
par(mfrow=c(3,4))
hist(Bad_rows$GM.CSF..14.,main = "GM-CSF",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$IL.1.alpha..22.,main = "IL-1a",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$IL.RA..29.,main = "IL-Ra",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$IL.2..35.,main = "IL-2",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$IL.4..39.,main = "IL-4",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$IL.6..43.,main = "IL-6",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$IL.8..46.,main = "IL-8",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$IL.10..52.,main = "IL-10",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$IL.12..57.,main = "IL-12",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$IL.18..61.,main = "IL-18",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$TNF.alpha..76.,main = "TNF-a",xlab="Cytokine Presence",col = "red")

par(mfrow=c(2,4))
hist(Bad_rows$MMP.1..18.,main = "MMP-1",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$MMP.2..21.,main = "MMP-2",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$MMP.3..27.,main = "MMP-3",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$MMP.7..30.,main = "MMP-7",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$MMP.9..33.,main = "MMP-9",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$MMP.10..37.,main = "MMP-10",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$MMP.12..44.,main = "MMP-12",xlab="Cytokine Presence",col = "red")
hist(Bad_rows$MMP.13..48.,main = "MMP-13",xlab="Cytokine Presence",col = "red")

####Heatmap

Lum_data$damage_score_24 <- NA
Lum_data$damage_score_24[which(Lum_data$Samples == "35-124-L")] <- 1 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "35-142-R")] <- 2 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "35-099-R")] <- 3 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "33-147-R")] <- 3 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "36-023-R")] <- 5 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "46-151-R")] <- 6 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "35-072-R")] <- 7 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "46-046-R")] <- 7 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "48-052-L")] <- 10 #previously 7
Lum_data$damage_score_24[which(Lum_data$Samples == "47-114-L")] <- 8 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "48-091-R")] <- 12 #previously 9
Lum_data$damage_score_24[which(Lum_data$Samples == "38-020-L")] <- 10 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "47-023-L")] <- 13 #previously 10
Lum_data$damage_score_24[which(Lum_data$Samples == "34-049-L")] <- 11 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "35-077-L")] <- 11 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "35-047-L")] <- 11 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "37-057-R")] <- 11 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "36-160-R")] <- 12 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "37-032-L")] <- 12 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "34-142-L")] <- 11 #previously 12
Lum_data$damage_score_24[which(Lum_data$Samples == "37-081-L")] <- 12 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "36-136-L")] <- 12 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "35-070-L")] <- 13 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "35-106-L")] <- 13 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "36-120-R")] <- 13 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "35-097-R")] <- 15 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "48-024-R")] <- 18 #previously 15
Lum_data$damage_score_24[which(Lum_data$Samples == "35-022-R")] <- 16 #good

Lum_data$damage_score_24[which(Lum_data$Samples == "45-146-R")] <- 6 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "48-064-L")] <- 13 #previously 10
Lum_data$damage_score_24[which(Lum_data$Samples == "38-055-R")] <- 5 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "47-151-L")] <- 6 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "48-017-L")] <- 5 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "34-156-R")] <- 15 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "35-068-L")] <- 11 #good
Lum_data$damage_score_24[which(Lum_data$Samples == "47-126-R")] <- 14 #good

library(pheatmap)
library(viridis)
library(scales)

### Heatmap with all samples, all MMPs/cytokines
# Specify colors
library(stringr)
colnames(Lum_data)[4:22] <- str_replace(unlist(strsplit(colnames(Lum_data)[4:22],"\\.\\."))[c(TRUE,FALSE)],"\\.","-")
colnames(Lum_data)[7] <- "IL-1α"
colnames(Lum_data)[9] <- "IL-1RA"
colnames(Lum_data)[22] <- "TNFα"
Lum_data$Cartilage_Outcome <- NA
Lum_data$Cartilage_Outcome[which(Lum_data$damage_score_24 >= 9)] <- "Bad"
#Lum_data$Cartilage_Outcome[which(Lum_data$damage_score_24 <= 14 & Lum_data$damage_score_24 >= 10)] <- "Fair"
Lum_data$Cartilage_Outcome[which(Lum_data$damage_score_24 <= 8)] <- "Good"
damage_colors <- viridis_pal(option="magma",begin = 1, end = 0.7)(16)
Lum_data$Treatment[Lum_data$Treatment=="RECON"] <- "ACLR"
Lum_data$Treatment[Lum_data$Treatment=="BE-REPAIR"] <- "BE-R"
annotation_colors = list(Treatment = c("ACLT"="mistyrose",
                                       "ACLR"="plum",
                                       "BE-R"="orchid4"),
                         Timepoint = c("Pre-op"="gray100",
                                       "1week"="gray90",
                                       "4weeks"="gray80",
                                       "12weeks"="gray70",
                                       "26weeks"="gray60",
                                       "52weeks"="gray50"),
                         damage_score_24 = c("1"=damage_colors[1],"2"=damage_colors[2],"3"=damage_colors[3],"4"=damage_colors[4],
                                             "5"=damage_colors[5],"6"=damage_colors[6],"7"=damage_colors[7],"8"=damage_colors[8],
                                             "9"=damage_colors[9],"10"=damage_colors[10],"11"=damage_colors[11],"12"=damage_colors[12],
                                             "13"=damage_colors[13],"14"=damage_colors[14],"15"=damage_colors[15],"16"=damage_colors[16]),
                         Cartilage_Outcome = c("Bad"="tomato","Good"="palegreen"))
pheatmap(t(Lum_data[,c(4:22)]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 16,
         fontsize_row = 18, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2")

#Heatmap with all samples, just MMPs
pheatmap(t(Lum_data[,c(grep("MMP",colnames(Lum_data)))]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2")
#Heatmaps with all samples just cytokines
pheatmap(t(Lum_data[,c(setdiff(c(4:22),grep("MMP",colnames(Lum_data))))]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2")
#Heatmap with Pre-op samples, all MMPs/cytokines
pheatmap(t(Lum_data[which(Lum_data$Timepoint=="Pre-op"),c(4:22)]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2")
#Heatmap with 1week samples, all MMPs/cytokines
pheatmap(t(Lum_data[which(Lum_data$Timepoint=="1week"),c(4:22)]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2")
#Heatmap with 4week samples, all MMPs/cytokines
pheatmap(t(Lum_data[which(Lum_data$Timepoint=="4weeks"),c(4:22)]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2")
#Heatmap with 12week samples, all MMPs/cytokines
pheatmap(t(Lum_data[which(Lum_data$Timepoint=="12weeks"),c(4:22)]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2")
#Heatmap with 26week samples, all MMPs/cytokines
pheatmap(t(Lum_data[which(Lum_data$Timepoint=="26weeks"),c(4:22)]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2")
#Heatmap with 52week samples, all MMPs/cytokines
pheatmap(t(Lum_data[which(Lum_data$Timepoint=="52weeks"),c(4:22)]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2")
#Heatmap with ACLT only
pheatmap(t(Lum_data[which(Lum_data$Treatment=="ACLT"),c(6,16,12,14)]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2", cluster_rows = FALSE)
#Heatmap with RECON only
pheatmap(t(Lum_data[which(Lum_data$Treatment=="RECON"),c(6,16,12,14)]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2", cluster_rows = FALSE)
#Heatmap with BE-R only
pheatmap(t(Lum_data[which(Lum_data$Treatment=="BE-REPAIR"),c(6,16,12,14)]), color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2", cluster_rows = FALSE)

pheatmap(t(Lum_data[which(Lum_data$Timepoint=="52weeks" | Lum_data$Timepoint=="4weeks" | Lum_data$Timepoint=="12weeks" | Lum_data$Timepoint=="26weeks"),c(4:22)]),
         color=c("mistyrose",damage_colors[16]), annotation_col=Lum_data[,c(2:3,29,30)], fontsize = 10,
         fontsize_row = 10, height=20, annotation_colors=annotation_colors,scale = "none", show_colnames=FALSE,legend_breaks=c(0,1),
         clustering_method = "ward.D2")

#Make dummy variables for Treatment groups
Lum_data$Repair <- 0
Lum_data$Repair[grep("BE-REPAIR",Lum_data$Treatment)] <- 1
Lum_data$Recon <- 0
Lum_data$Recon[grep("RECON",Lum_data$Treatment)] <- 1
Lum_data$ACLT <- 0
Lum_data$ACLT[grep("ACLT",Lum_data$Treatment)] <- 1

#Adding Gait data
setwd("") # <---- insert your own path to the gait data file here
Gait_raw <- read_excel("gait_raw_data_for_PLOS_ONE.xlsx")
#Isolating Gait data from Surgical limbs, no Gait data collected for 1 wk
Gait_contra <- Gait_raw[which(Gait_raw$`Code-Limb` == "CONTRALATERAL"),]
Gait_raw <- Gait_raw[which(Gait_raw$`Code-Limb` == "SURGICAL"),]
table(Gait_raw$Timepoint)
table(Gait_contra$Timepoint)
table(Lum_data$Timepoint)
Gait_4_week <- Gait_raw[which(Gait_raw$Timepoint == "4weeks"),]
Gait_12_week <- Gait_raw[which(Gait_raw$Timepoint == "12weeks"),]
Gait_26_week <- Gait_raw[which(Gait_raw$Timepoint == "26weeks"),]
Gait_52_week <- Gait_raw[which(Gait_raw$Timepoint == "52weeks"),]
Gait_4_week_contra <- Gait_contra[which(Gait_contra$Timepoint == "4weeks"),]
Gait_12_week_contra <- Gait_contra[which(Gait_contra$Timepoint == "12weeks"),]
Gait_26_week_contra <- Gait_contra[which(Gait_contra$Timepoint == "26weeks"),]
Gait_52_week_contra <- Gait_contra[which(Gait_contra$Timepoint == "52weeks"),]
#Put Gait_contra rows in same order as Gait_raw
Gait_4_week_contra <- Gait_4_week_contra[order(match(Gait_4_week$ID,Gait_4_week_contra$ID)),]
Gait_12_week_contra <- Gait_12_week_contra[order(match(Gait_12_week$ID,Gait_12_week_contra$ID)),]
Gait_26_week_contra <- Gait_26_week_contra[order(match(Gait_26_week$ID,Gait_26_week_contra$ID)),]
Gait_52_week_contra <- Gait_52_week_contra[order(match(Gait_52_week$ID,Gait_52_week_contra$ID)),]
Lum_data$stride_length_hind <- NA
Lum_data$stride_length_contra <- NA
Lum_data$stride_length_ratio <- NA
Lum_data$stride_vel_hind <- NA
Lum_data$stride_vel_contra <- NA
Lum_data$stride_vel_ratio <- NA
Lum_data$Max_pressure_hind <- NA
Lum_data$Max_pressure_contra <- NA
Lum_data$Max_pressure_ratio <- NA
Lum_data$Max_force_BW_hind <- NA
Lum_data$Max_force_BW_contra <- NA
Lum_data$Max_force_BW_ratio <- NA
Lum_data$Impulse_BW_hind <- NA
Lum_data$Impulse_BW_contra <- NA
Lum_data$Impulse_BW_ratio <- NA
for(i in 1:dim(Gait_4_week)[1]){
  Lum_data$stride_length_hind[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week$`ID-Limb`[i])] <- Gait_4_week$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_hind[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week$`ID-Limb`[i])] <- Gait_4_week$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_hind[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week$`ID-Limb`[i])] <- Gait_4_week$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_hind[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week$`ID-Limb`[i])] <- Gait_4_week$`Maximum Force (%BW) Hind`[i]
  Lum_data$Impulse_BW_hind[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week$`ID-Limb`[i])] <- Gait_4_week$`Impulse (%BW*sec) Hind`[i]
  
  Lum_data$stride_length_contra[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week_contra$`ID-Limb`[i])] <- Gait_4_week_contra$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_contra[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week_contra$`ID-Limb`[i])] <- Gait_4_week_contra$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_contra[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week_contra$`ID-Limb`[i])] <- Gait_4_week_contra$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_contra[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week_contra$`ID-Limb`[i])] <- Gait_4_week_contra$`Maximum Force (%BW) Hind`[i]
  Lum_data$Impulse_BW_contra[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week$`ID-Limb`[i])] <- Gait_4_week_contra$`Impulse (%BW*sec) Hind`[i]
  
  Lum_data$stride_length_ratio[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week$`ID-Limb`[i])] <- Gait_4_week$`Stride Length (cm) Hind`[i]/Gait_4_week_contra$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_ratio[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week$`ID-Limb`[i])] <- Gait_4_week$`Stride Velocity (cm/sec) Hind`[i]/Gait_4_week_contra$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_ratio[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week$`ID-Limb`[i])] <- Gait_4_week$`Maximum Peak Pressure (KPa) Hind`[i]/Gait_4_week_contra$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_ratio[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week$`ID-Limb`[i])] <- as.numeric(Gait_4_week$`Maximum Force (%BW) Hind`[i])/as.numeric(Gait_4_week_contra$`Maximum Force (%BW) Hind`[i])
  Lum_data$Impulse_BW_ratio[which(Lum_data$Timepoint == "4weeks" & Lum_data$Samples == Gait_4_week$`ID-Limb`[i])] <- as.numeric(Gait_4_week$`Impulse (%BW*sec) Hind`[i])/as.numeric(Gait_4_week_contra$`Impulse (%BW*sec) Hind`[i])
}
for(i in 1:dim(Gait_12_week)[1]){
  Lum_data$stride_length_hind[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week$`ID-Limb`[i])] <- Gait_12_week$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_hind[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week$`ID-Limb`[i])] <- Gait_12_week$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_hind[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week$`ID-Limb`[i])] <- Gait_12_week$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_hind[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week$`ID-Limb`[i])] <- Gait_12_week$`Maximum Force (%BW) Hind`[i]
  Lum_data$Impulse_BW_hind[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week$`ID-Limb`[i])] <- Gait_12_week$`Impulse (%BW*sec) Hind`[i]
  
  Lum_data$stride_length_contra[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week_contra$`ID-Limb`[i])] <- Gait_12_week_contra$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_contra[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week_contra$`ID-Limb`[i])] <- Gait_12_week_contra$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_contra[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week_contra$`ID-Limb`[i])] <- Gait_12_week_contra$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_contra[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week_contra$`ID-Limb`[i])] <- Gait_12_week_contra$`Maximum Force (%BW) Hind`[i]
  Lum_data$Impulse_BW_contra[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week$`ID-Limb`[i])] <- Gait_12_week_contra$`Impulse (%BW*sec) Hind`[i]
  
  Lum_data$stride_length_ratio[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week$`ID-Limb`[i])] <- Gait_12_week$`Stride Length (cm) Hind`[i]/Gait_12_week_contra$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_ratio[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week$`ID-Limb`[i])] <- Gait_12_week$`Stride Velocity (cm/sec) Hind`[i]/Gait_12_week_contra$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_ratio[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week$`ID-Limb`[i])] <- Gait_12_week$`Maximum Peak Pressure (KPa) Hind`[i]/Gait_12_week_contra$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_ratio[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week$`ID-Limb`[i])] <- as.numeric(Gait_12_week$`Maximum Force (%BW) Hind`[i])/as.numeric(Gait_12_week_contra$`Maximum Force (%BW) Hind`[i])
  Lum_data$Impulse_BW_ratio[which(Lum_data$Timepoint == "12weeks" & Lum_data$Samples == Gait_12_week$`ID-Limb`[i])] <- as.numeric(Gait_12_week$`Impulse (%BW*sec) Hind`[i])/as.numeric(Gait_12_week_contra$`Impulse (%BW*sec) Hind`[i])
}
for(i in 1:dim(Gait_26_week)[1]){
  Lum_data$stride_length_hind[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week$`ID-Limb`[i])] <- Gait_26_week$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_hind[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week$`ID-Limb`[i])] <- Gait_26_week$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_hind[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week$`ID-Limb`[i])] <- Gait_26_week$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_hind[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week$`ID-Limb`[i])] <- Gait_26_week$`Maximum Force (%BW) Hind`[i]
  Lum_data$Impulse_BW_hind[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week$`ID-Limb`[i])] <- Gait_26_week$`Impulse (%BW*sec) Hind`[i]
  
  Lum_data$stride_length_contra[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week_contra$`ID-Limb`[i])] <- Gait_26_week_contra$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_contra[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week_contra$`ID-Limb`[i])] <- Gait_26_week_contra$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_contra[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week_contra$`ID-Limb`[i])] <- Gait_26_week_contra$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_contra[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week_contra$`ID-Limb`[i])] <- Gait_26_week_contra$`Maximum Force (%BW) Hind`[i]
  Lum_data$Impulse_BW_contra[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week$`ID-Limb`[i])] <- Gait_26_week_contra$`Impulse (%BW*sec) Hind`[i]
  
  Lum_data$stride_length_ratio[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week$`ID-Limb`[i])] <- Gait_26_week$`Stride Length (cm) Hind`[i]/Gait_26_week_contra$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_ratio[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week$`ID-Limb`[i])] <- Gait_26_week$`Stride Velocity (cm/sec) Hind`[i]/Gait_26_week_contra$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_ratio[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week$`ID-Limb`[i])] <- Gait_26_week$`Maximum Peak Pressure (KPa) Hind`[i]/Gait_26_week_contra$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_ratio[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week$`ID-Limb`[i])] <- as.numeric(Gait_26_week$`Maximum Force (%BW) Hind`[i])/as.numeric(Gait_26_week_contra$`Maximum Force (%BW) Hind`[i])
  Lum_data$Impulse_BW_ratio[which(Lum_data$Timepoint == "26weeks" & Lum_data$Samples == Gait_26_week$`ID-Limb`[i])] <- as.numeric(Gait_26_week$`Impulse (%BW*sec) Hind`[i])/as.numeric(Gait_26_week_contra$`Impulse (%BW*sec) Hind`[i])
}
for(i in 1:dim(Gait_52_week)[1]){
  Lum_data$stride_length_hind[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week$`ID-Limb`[i])] <- Gait_52_week$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_hind[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week$`ID-Limb`[i])] <- Gait_52_week$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_hind[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week$`ID-Limb`[i])] <- Gait_52_week$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_hind[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week$`ID-Limb`[i])] <- Gait_52_week$`Maximum Force (%BW) Hind`[i]
  Lum_data$Impulse_BW_hind[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week$`ID-Limb`[i])] <- Gait_52_week$`Impulse (%BW*sec) Hind`[i]
  
  Lum_data$stride_length_contra[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week_contra$`ID-Limb`[i])] <- Gait_52_week_contra$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_contra[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week_contra$`ID-Limb`[i])] <- Gait_52_week_contra$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_contra[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week_contra$`ID-Limb`[i])] <- Gait_52_week_contra$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_contra[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week_contra$`ID-Limb`[i])] <- Gait_52_week_contra$`Maximum Force (%BW) Hind`[i]
  Lum_data$Impulse_BW_contra[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week$`ID-Limb`[i])] <- Gait_52_week_contra$`Impulse (%BW*sec) Hind`[i]
  
  Lum_data$stride_length_ratio[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week$`ID-Limb`[i])] <- Gait_52_week$`Stride Length (cm) Hind`[i]/Gait_52_week_contra$`Stride Length (cm) Hind`[i]
  Lum_data$stride_vel_ratio[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week$`ID-Limb`[i])] <- Gait_52_week$`Stride Velocity (cm/sec) Hind`[i]/Gait_52_week_contra$`Stride Velocity (cm/sec) Hind`[i]
  Lum_data$Max_pressure_ratio[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week$`ID-Limb`[i])] <- Gait_52_week$`Maximum Peak Pressure (KPa) Hind`[i]/Gait_52_week_contra$`Maximum Peak Pressure (KPa) Hind`[i]
  Lum_data$Max_force_BW_ratio[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week$`ID-Limb`[i])] <- as.numeric(Gait_52_week$`Maximum Force (%BW) Hind`[i])/as.numeric(Gait_52_week_contra$`Maximum Force (%BW) Hind`[i])
  Lum_data$Impulse_BW_ratio[which(Lum_data$Timepoint == "52weeks" & Lum_data$Samples == Gait_52_week$`ID-Limb`[i])] <- as.numeric(Gait_52_week$`Impulse (%BW*sec) Hind`[i])/as.numeric(Gait_52_week_contra$`Impulse (%BW*sec) Hind`[i])
}

#Removing pre-op data
Lum_data <- Lum_data[which(Lum_data$Timepoint != "Pre-op"),]
Lum_data <- Lum_data[which(Lum_data$Timepoint != "1week"),]

colnames(Lum_data)
for (i in 34:48) {
  print(colnames(Lum_data)[i])
  print(cor.test(Lum_data$ACLT,as.numeric(Lum_data[,i]),method = "kendall")[3])
  print(cor.test(Lum_data$Recon,as.numeric(Lum_data[,i]),method = "kendall")[3])
  print(cor.test(Lum_data$Repair,as.numeric(Lum_data[,i]),method = "kendall")[3])
}

for (i in 34:48) {
  Lum_data[,i] <- as.numeric(Lum_data[,i])
  Lum_data[,i] <- (Lum_data[,i]-mean(Lum_data[,i]))/(sd(c(Lum_data[,i])))
}

# example and usage of gaussian mixture models for clustering https://tinyheero.github.io/2015/10/13/mixture-model.html
#See if there is a separation in the damage distributions
qplot(Lum_data$damage_score_24, geom="histogram",bins=18,col=I("red"),main = "Histogram for Damage")
#Split into two gaussians
plot_mix_comps <- function(x, mu, sigma, lam) {
  lam * dnorm(x, mu, sigma)
}
mixmdl <- normalmixEM(Lum_data$damage_score_24, k = 2)
data.frame(x = mixmdl$x) %>%
  ggplot() +
  geom_histogram(aes(x, ..density..), binwidth = 1, colour = "black", 
                 fill = "white") +
  stat_function(geom = "line", fun = plot_mix_comps,
                args = list(mixmdl$mu[1], mixmdl$sigma[1], lam = mixmdl$lambda[1]),
                colour = "red", lwd = 1.5) +
  stat_function(geom = "line", fun = plot_mix_comps,
                args = list(mixmdl$mu[2], mixmdl$sigma[2], lam = mixmdl$lambda[2]),
                colour = "blue", lwd = 1.5) +
  ylab("Density") + xlab("Macroscopic Cartilage Damage Score") + theme(axis.title = element_text(size = 24)) +
  theme(axis.text = element_text(size = 22))   
post.df <- as.data.frame(cbind(x = mixmdl$x, mixmdl$posterior))
head(post.df,30)
#If the probability of a point belonging to group 1 (Good OA group) is > 0.5, then it goes to group 1.
post.df %>%
  mutate(label = ifelse(comp.1 > 0.5, 1, 2)) %>% 
  ggplot(aes(x = factor(label))) +
  geom_bar() +
  xlab("Component") +
  ylab("Number of Data Points")
rownames(Lum_data) <- rownames(post.df)
Lum_data$Binary_OA_outcomes <- NA
Lum_data$Binary_OA_outcomes[which(post.df$comp.1 > 0.5)] <- 1
Lum_data$Binary_OA_outcomes[which(post.df$comp.2 >= 0.5)] <- 0
#The split is [1-8] = Good, [9-16] = Bad
table(Lum_data$Binary_OA_outcomes[random_rows[1:91]])

set.seed(42)
random_rows <- sample(nrow(Lum_data))
colnames(Lum_data)
#30% going into cold storage for testing
setwd("") # <---- choose your own directory data storage
#Training sets
write.csv(Lum_data[random_rows[1:91],c(2,4:21,34:48,49)],"Luminex_data_training_set_for_LOOCV.csv") #SF + gait
write.csv(Lum_data[random_rows[1:91],c(2,4:21,49)],"Luminex_data_training_set_for_LOOCV.csv") #SF
write.csv(Lum_data[random_rows[1:91],c(2,34:48,49)],"Luminex_data_training_set_for_LOOCV.csv") #gait
write.csv(Lum_data[random_rows[1:91],c(2,7,9,12,14,15,19,6,8,16,18,40,41,49)],"Luminex_data_training_set_for_LOOCV.csv") #GEE
#26 subjects (92 profiles) in training
#sample_split <- unique(sample(36))
training_set <- Lum_data[which(Lum_data$Samples %in% unique(Lum_data$Samples)[sample_split[1:26]]),]
write.csv(training_set[,c(2,4:21,34:48,49)],"Luminex_data_training_set_for_LOOCV.csv") #SF + gait
write.csv(training_set[,c(2,4:21,49)],"Luminex_data_training_set_for_LOOCV.csv") #SF
write.csv(training_set[,c(2,7,9,12,14,15,19,6,8,16,18,40,41,49)],"Luminex_data_training_set_for_LOOCV.csv") #GEE
#Testing sets
write.csv(Lum_data[random_rows[92:129],c(2,4:21,34:48,49)],"Luminex_data_testing_set_for_LOOCV.csv") #SF + gait 33 C = 10000.0, 0.1
write.csv(Lum_data[random_rows[92:129],c(2,4:21,49)],"Luminex_data_testing_set_for_LOOCV.csv") #SF 18 C = 1000.0
write.csv(Lum_data[random_rows[92:129],c(2,34:48,49)],"Luminex_data_testing_set_for_LOOCV.csv") #gait 15 C = 0.1
write.csv(Lum_data[random_rows[92:129],c(2,7,9,12,14,15,19,6,8,16,18,40,41,49)],"Luminex_data_testing_set_for_LOOCV.csv") #GEE 12 1.0
#10 subjects (37 profiles) in testing
testing_set <- Lum_data[which(Lum_data$Samples %in% unique(Lum_data$Samples)[sample_split[27:36]]),]
write.csv(testing_set[,c(2,4:21,34:48,49)],"Luminex_data_testing_set_for_LOOCV.csv") #SF + gait 33 C = 10000.0, 0.1
write.csv(testing_set[,c(2,4:21,49)],"Luminex_data_testing_set_for_LOOCV.csv") #SF 18 C = 1000.0
write.csv(testing_set[,c(2,7,9,12,14,15,19,6,8,16,18,40,41,49)],"Luminex_data_testing_set_for_LOOCV.csv") #GEE 12 1.0

summary(glm(Binary_OA_outcomes ~ GM.CSF..14.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ TNF.alpha..76.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ IL.1.alpha..22.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ IL.RA..29.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ IL.2..35.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ IL.4..39.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ IL.6..43.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ IL.8..46.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ IL.10..52.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ IL.12..57.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ IL.18..61.,family = binomial(link="logit"),data = Lum_data))[12]

summary(glm(Binary_OA_outcomes ~ MMP.1..18.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ MMP.2..21.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ MMP.3..27.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ MMP.7..30.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ MMP.9..33.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ MMP.10..37.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ MMP.12..44.,family = binomial(link="logit"),data = Lum_data))[12]
summary(glm(Binary_OA_outcomes ~ MMP.13..48.,family = binomial(link="logit"),data = Lum_data))[12]

summary(glmer(Binary_OA_outcomes ~ GM.CSF..14. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ TNF.alpha..76. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.1.alpha..22. + (1|Samples), data=Lum_data, family=binomial(link="logit")))
summary(glmer(Binary_OA_outcomes ~ IL.RA..29. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.2..35. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.4..39. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.6..43. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.8..46. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.10..52. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.12..57. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.18..61. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]

summary(glmer(Binary_OA_outcomes ~ MMP.1..18. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.2..21. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.3..27. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.7..30. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.9..33. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.10..37. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.12..44. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.13..48. + (1|Samples), data=Lum_data, family=binomial(link="logit")))[10]

summary(glmer(Binary_OA_outcomes ~ GM.CSF..14. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ TNF.alpha..76. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.1.alpha..22. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.RA..29. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.2..35. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.4..39. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.6..43. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.8..46. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.10..52. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.12..57. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ IL.18..61. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]

summary(glmer(Binary_OA_outcomes ~ MMP.1..18. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.2..21. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.3..27. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.7..30. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.9..33. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.10..37. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.12..44. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]
summary(glmer(Binary_OA_outcomes ~ MMP.13..48. + (1|Timepoint), data=Lum_data, family=binomial(link="logit")))[10]

Lum_data_gee <- Lum_data
Lum_data_gee <- Lum_data_gee[order(Lum_data_gee$Samples),]
Lum_data_gee$Samples <- factor(Lum_data_gee$Samples)
(pnorm(summary(gee(Binary_OA_outcomes ~ GM.CSF..14., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(pnorm(summary(gee(Binary_OA_outcomes ~ TNF.alpha..76., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ IL.1.alpha..22., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ IL.RA..29., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ IL.2..35., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ IL.4..39., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ IL.6..43., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ IL.8..46., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ IL.10..52., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ IL.12..57., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ IL.18..61., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2

(1-pnorm(summary(gee(Binary_OA_outcomes ~ MMP.1..18., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ MMP.2..21., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ MMP.3..27., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ MMP.7..30., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ MMP.9..33., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ MMP.10..37., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ MMP.12..44., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ MMP.13..48., id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2

(1-pnorm(summary(gee(Binary_OA_outcomes ~ stride_length_hind, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ stride_length_contra, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(pnorm(summary(gee(Binary_OA_outcomes ~ stride_length_ratio, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ stride_vel_hind, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ stride_vel_contra, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(pnorm(summary(gee(Binary_OA_outcomes ~ stride_vel_ratio, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(pnorm(summary(gee(Binary_OA_outcomes ~ Max_pressure_hind, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(pnorm(summary(gee(Binary_OA_outcomes ~ Max_pressure_contra, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(pnorm(summary(gee(Binary_OA_outcomes ~ Max_pressure_ratio, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ Max_force_BW_hind, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ Max_force_BW_contra, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(pnorm(summary(gee(Binary_OA_outcomes ~ Max_force_BW_ratio, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(pnorm(summary(gee(Binary_OA_outcomes ~ Impulse_BW_hind, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(pnorm(summary(gee(Binary_OA_outcomes ~ Impulse_BW_contra, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2
(1-pnorm(summary(gee(Binary_OA_outcomes ~ Impulse_BW_ratio, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[6]))*2

coef <- summary(gee(Binary_OA_outcomes ~ Impulse_BW_ratio, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[2]
se <- summary(gee(Binary_OA_outcomes ~ Impulse_BW_ratio, id=Samples, data=Lum_data_gee, family=binomial))$coefficients[4]
paste0("coef: ",signif(coef, digits = 4)," OR: ",signif(exp(coef),digits = 4))
paste0("95% CI: ",signif(exp(c(coef-(1.96*se),coef+(1.96*se))),digits = 4))